//
// Distributions.cxx
//
// Created by TB on 11/30/11.
// Copyright 2011 BNL. All rights reserved.
//

#include <iostream>

#include <boost/foreach.hpp>
#include <boost/shared_ptr.hpp>

#include "tppmc/CollinsFragmentation.h"
#include "tppmc/Configuration.h"
#include "tppmc/ConstantFragmentation.h"
#include "tppmc/Distributions.h"
#include "tppmc/DssFragmentation.h"
#include "tppmc/Dssv.h"
#include "tppmc/CollinsKt.h"
#include "tppmc/Grv98Pdf.h"
#include "tppmc/sivers/Torino2009.h"
#include "tppmc/sivers/Positivity.h"
#include "tppmc/SofferBound.h"
#include "tppmc/SofferBoundProkudin.h"
#include "tppmc/TransversityProkudin.h"

#ifdef USING_LHAPDF
#include "tppmc/LhapdfDistribution.h"
#endif

namespace tppmc {

   typedef boost::shared_ptr<PartonDistribution> PartonPtr;
   typedef boost::shared_ptr<FragmentationFunction> FragPtr;

   Distributions::Distributions() {
      // Assign PDF set. No option currently.
      mUnpolarised = PartonPtr(new Grv98Pdf(Grv98Pdf::LO));
      // Assign Sivers set. Options are Anselmino extraction or
      // positivity bound (needs to be implemented).
      mSivers = PartonPtr(new sivers::Positivity());
      // Assign transversity set. Options are Prokudin extraction
      // or Soffer bound.
      mTransversity =
      PartonPtr(new SofferBoundProkudin(SofferBoundProkudin::LO));
      // Assign fragmentation function. Options are constant fraction
      // of parton momentum, DSS, DSS + Collins, KKP and AKK (KKP
      // and AKK need to be implemented).
      mFragmentation = FragPtr(new ConstantFragmentation(1.));
   }

   Distributions::Distributions(const Configuration& config) {
#ifdef USING_LHAPDF
      if(config.LhapdfSet() > 0) {
         mUnpolarised = PartonPtr(new LhapdfDistribution(
            config.LhapdfSet(), config.LhapdfMember())
         );
      } // if
      else {
         mUnpolarised = PartonPtr(new Grv98Pdf(Grv98Pdf::LO));
      } // else
#else
      mUnpolarised = PartonPtr(new Grv98Pdf(Grv98Pdf::LO));
#endif
      // Assign Sivers set. Only choice is Anselmino extraction or-
      // positivity bound based on the unpolarised distribution.
      switch(config.SiversType()) {
         case 1: // Anselmino et al. 2009, normal (not positivity bound)
            mSivers = PartonPtr(new sivers::Torino2009(false));
            break;
         case 2: // Positivity bound based on unpolarised PDF
            mSivers = PartonPtr(new sivers::Positivity(mUnpolarised));
            break;
         default: // No sivers
            mSivers = PartonPtr(new NullDistribution);
            break;
      } // switch

      // Assign transversity set. Options are Prokudin extraction
      // or Soffer bound.
      /** \todo Add command line options to switch between transversity sets,
      select sea quarks 0 or Soffer bound etc*/
      PartonPtr dssv(new Dssv);
      if(config.UseTransversitySofferBound()) {
         mTransversity = PartonPtr(
            new transversity::SofferBound(mUnpolarised, dssv, false)
         );
      } // if
      else {
         mTransversity = PartonPtr(new NullDistribution);
      } // else
      
      // Assign fragmentation function.
      // If fragmentation is not requested hadrons are generated with a
      // constant fraction of the parton momentum.
      // Options for fragmentation are DSS, DSS + Collins, KKP and AKK (KKP
      // and AKK need to be implemented).
      switch(config.FragmentationFunction()) {
         case 0: // No fragmentation
            mFragmentation = FragPtr(new NullFragmentation());
            break;
         case 1: // DSS
            mFragmentation = FragPtr(
               new DssFragmentation(config.Order(), config.MeanJt2())
            );
            break;
         case 2: // Constant fragmentation
            mFragmentation = FragPtr(new ConstantFragmentation(1.));
            break;
         case 4: // DSS + Collins
            mFragmentation = FragPtr(new CollinsFragmentation);
            (dynamic_cast<CollinsFragmentation*>(mFragmentation.get()))
            ->SetPositivityBound(config.UseCollinsPositivityBound());
            break;
         case 5: // DSS+Collins with kT-dependent positivity bound
            mFragmentation = FragPtr(new collins::PositivityKt(
               boost::shared_ptr<FragmentationFunction>(
               new DssFragmentation(config.Order(), config.MeanJt2())),
               config.MeanJt2()));
            break;
         default:
            mFragmentation = FragPtr(new ConstantFragmentation(1.));
            break;
      } // switch
      mUnpolarised->Initialise();
      mSivers->Initialise();
      mTransversity->Initialise();
      mFragmentation->Initialise();
   }

   Distributions::~Distributions() {
   }

   PartonPtr Distributions::Unpolarised() const {
      return mUnpolarised;
   }
   
   PartonPtr Distributions::Sivers() const {
      return mSivers;
   }
   
   PartonPtr Distributions::Transversity() const {
      return mTransversity;
   }
   
   FragPtr Distributions::Fragmentation() const {
      return mFragmentation;
   }
   
   std::list<PartonPtr> Distributions::Pdfs() const {
      std::list<PartonPtr> pdfs;
      
      if(mUnpolarised) pdfs.push_back(mUnpolarised);
      if(mSivers) pdfs.push_back(mSivers);
      if(mTransversity) pdfs.push_back(mTransversity);
      
      return pdfs;
   }

   void Distributions::Speak() const {
      std::cout << "Unpolarised:   " << Unpolarised()->Name()   << std::endl;
      std::cout << "Sivers:        " << Sivers()->Name()        << std::endl;
      std::cout << "Transversity:  " << Transversity()->Name()  << std::endl;
      std::cout << "Fragmentation: " << Fragmentation()->Name() << std::endl;
   }
   
} // namespace tppmc
